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Abstract 



This paper deals with the parallel simulation of delamination problems at the meso- 
scale by means of multi-scale methods, the aim being the Virtual Delamination Testing 
of Composite parts. In the non-linear context, Domain Decomposition Methods are 
mainly used as a solver for the tangent problem to be solved at each iteration of a 
Newton-Raphson algorithm. In case of strongly non linear and heterogeneous problems, 

^^ this procedure may lead to severe difficulties. The paper focuses on methods to circum- 

^ vent these problems, which can now be expressed using a relatively general framework, 

even though the different ingredients of the strategy have emerged separately. We rely 
here on the micro- macro framework proposed in [36J. The method proposed in this 
paper introduces three additional features: (i) the adaptation of the macro-basis to sit- 
uations where classical homogenization does not provide a good preconditioner, (ii) the 
use of non-linear relocalization to decrease the number of global problems to be solved 
in the case of unevenly distributed non-linearities, (iii) the adaptation of the approxi- 
mation of the local Schur complement which governs the convergence of the proposed 
IV iterative technique. Computations of delamination and delamination-buckling interac- 

. ^h tion with contact on potentially large delaminated areas are used to illustrate those 

^ aspects. 



1 Introduction 

Despite the many studies since the beginning of the eighties, the prediction of delami- 
nation remains a challenge, both from a scientific and from an industrial point of view. 
Many difficulties encountered in this field are inherent to composite modeling. The in- 
terested reader will find a survey of the literature devoted to the modeling of composites 
in |26j. Simulating delamination requires to tackle specific difficulties (i) the complex 



state of stress which leads to the initiation and the propagation of delamination (ii) 
the size and the complexity of the problem to be solved when dealing with composite 
structures. One of the industrial objectives is to replace some of the numerous physical 
tests required to assess the damage tolerance of composites by numerical simulations 
|4J, or in other words to perform Virtual Delamination Testing (VDT). 

It is today well accepted at an industrial level that a realistic VDT procedure should 
rely on a meso-scale description of the laminates. Such a fine-scale description requires 
a Finite Element discretisation whose characteristic size is at most the thickness of 
individual plies, which is of the order of a tenth of a millimetre. For the treatment 
of industrial parts, the size of the resulting model is therefore huge, which leads to 
considering the applicability of parallel computations to those large and complex non- 
linear cases. Indeed, these problems potentially involve multiple delamination fronts, 
buckling or contact, and therefore instabilities. 

The paper describes some attempts to deal with such problems within a multi-scale 
parallel framework when delamination is the only damage mode that is taken into 
account. Delamination is here modeled by means of a cohesive zone model. Several 
such models have been proposed in the literature, for example in [6, 52, 49, 46] [TO]. 
In this paper we make use of the model proposed in [5 J but similar results could be 
obtained with other models once calibrated. The main features of the model are briefly 
presented in Section [2] 

In the remainder of the paper, the basics of the multi-scale strategy are described 
and test cases are presented to demonstrate its numerical efficiency. The contribution 
of the paper is the definition of a framework that encompasses different aspects of the 
strategy that have emerged separately in [341 |25j [13j |29l [3l [47]. Indeed, our previous 
publications focussed on particular issues related to the present topic. It now appears 
that the set of tools that we have developed in the past can help formulate a general 
non-linear multi-scale approach that is adapted to the treatment of large laminated and 
damageable structures taking into account possible buckling and contact. Multi-scale 
refers here to a computational technique which involves different scales and which aims 
at finding the solution over the whole structure at the smallest scale of interest. In this 
domain, one of the main issues to be addressed is the transfer of information from one 
scale to another. 

A large number of multiscale methods rely on classical homogenisation of heteroge- 
neous media, which was initiated by [27J. For linear periodic media the most efficient 
method was initiated in [48]. Further developments and related computational ap- 
proaches can be found in [53l [201 SSI EU [19] . An overview of these methods can 
be found in [54J. In computational homogenization approaches the resolution of the 
macro-problem leads to effective values of the unknowns; then, the micro solution is 
calculated locally based on the macro solution. The fundamental assumption, besides 
periodicity, is that the ratio of the characteristic length of the small scale to the char- 
acteristic length of the large scale must be small. All areas which do not satisfy the 
hypothesis of scale separation such as boundary zones or vicinity of cracks, where the 
material cannot be homogenized, require special treatment. Therefore a large amount 
of recent studies aim at applying enhanced homogenization schemes to localized failure, 
see for instance [i^ l7|[T^ ] . These extension, which are dedicated to overcome some of 
the difficulties for a given type of problem are not suited yet to deal with edge effects 



and large cracks and therefore not directly suitable to the treatment of delamination. 

The proposed method is based on a Domain Decomposition approach. These meth- 
ods have been developed initially for large linear problems for which they can now be 
considered mature. The most well-known methods are FETI [18J and BDD [41J which 
where designed for perfect interface and extended, in the case of FETI, to contact (with 
or without friction) [16] or prescribed displacement gaps [32J. FETI and BDD rely on 
a static condensation for each sub-domain. The resulting interface problem, which is 
large, is solved iteratively using a preconditioner which is associated to the balance of 
the rigid-body momentum of the sub-domains. In the case of heterogeneous media, the 
macro-space of rigid body is not necessarily the most appropriate one and a specific 
deflation can be performed [17] . It is therefore interesting to define the macro-space as 
the space which is associated with the homogenized response of the subdomains. The 
use of homogenization in the context of a Domain Decomposition Method has been 
proposed in the micro-macro approach of [35J. 

Homogenization being a natural tool for composites, we base our approach on this 
micro-macro framework. The decomposition of the structure into subdomains is chosen 
to be compatible with the mesomodeling of the laminate. The micro-macro method 
relies on mixed interface conditions between subdomains, which allows to treat any 
kind of interfacial constitutive law directly at the interface level [34] , Three types 
of interfaces are considered here: damageable interfaces and contact interfaces that 
are compatible with the mesomodeling, and perfect interfaces within the plies to split 
the problem into smaller pieces that are suited to parallelism. This leads to a large 
number of substructures, and therefore to a large macro-problem. A third scale is 
then introduced to perform a low cost and precise approximation of the homogeneized 
macro-problem by using a BDD preconditioner [29J, as described in Section [3] Note 
that in this paper only unilateral contact is considered but the framework is well suited 
to friction [11] . 

In the context of non-linear problems, Domain Decomposition Methods are mainly 
used as a solver for the tangent problem to be solved at each iteration of a Newton- 
Raphson based algorithm [38] |22] |30] , with if required an arc-length control. Several 
problems can be encountered in this case. One of such problems is that the choice of 
the macro-space, and the implied homogenized behavior, is in general not adapted to a 
situation where a crack interacts with a subdomain. In other words, a preconditioner 
designed for a problem without crack can become ineffective in the presence of a crack. 
It has been shown in [25, 24J, that for cracks crossing interfaces of the domain decom- 
position, a modification of the macro-basis that adds a few degrees of freedom only 
makes it possible to obtain the same convergence rate for a cracked domain and for an 
un-cracked one. Unfortunately such an adaptation is not possible in the case of cracks 
propagating at the interface between subdomains. We have therefore followed another 
line, the one of non-linear relocalization [13, 44J. This approach aims at solving the 
problem of the strong deterioration of the convergence rate due to localized non-linear 
effect [9J. Following this line of thoughts, we have applied a non-linear relocalization 
approach on the set of subdomains that are connected to the front of delamination, 
which is described in Section [4][3j. Note that our approach naturally avoids to conduct 
useless nonlinear computations far from the delamination front which is the goal of the 
prediction procedure developped in [39J for Newton- Schur strategies. 



Having in mind the problem of the tolerance of laminates to compression after im- 
pact [TS] we have started to work on the question of the interaction between potentially 
large growing delaminated areas and buckling, taking into account the contact condi- 
tion [47]. In order to reconnect the micro part of the solution, the iterative technique 
introduces parameters which are associated to the influence of neighboring subdomains. 
These so called "search directions" correspond to Robin parameters whose values were 
optimized for 3D structures. These parameters must be adapted to the aspect ratios of 
the slender structures by introducing well-chosen anisotropic coefficients. Moreover, in 
the case of contact between slender structures over large delaminated areas, the search 
directions should be optimized according to the interface's state (open or closed) in 
order to obtain a reasonable convergence rate of the iterative solver. These points are 
discussed in Section |5j However, changing these parameters requires a refactorization 
of the homogenized operators, which is computationally expensive. A compromise be- 
tween convergence speed and unit cost of iterations must be found so that the total 
CPU time is minimised. 



2 The interface model 

The interface model relates the jump of displacements [u] = u , — u p to the normal 
Cauchy stress t — a .n = a r n on the interface V between two plies p and p' '. Here, n 
denotes the outer normal to ply p on T pp f . In this problem, large displacements are con- 
sidered but the jump of displacement is assumed to be small prior to full delamination. 
The tractions t are uniquely defined and related to the displacement discontinuities by 
means of the following damageable constitutive relation: 



t = K.[u] 



(1) 



The expression of the stiffness operator K^ can be given in the reference frame (iV 1 , N_ 2 •> n) 
(as depicted in Figure with n = N^p = "Ply 1" and p' = "Ply 2"): 



(2) 



< > + is here the positive indicator function. 

The local damage variables &i are introduced into the interface model in order to 
simulate its softening behavior when the structure is loaded. Their values range from 
(healthy interface point) to 1 (completely damaged interface point). The parameters 
d{ are related to the local energy release rates Y\ of the interface's degradation modes 
as follows: 



/(l- 


- d\)k\ 







\ 







(1 - d 2 )k 2 







V 








('■ 


- ^ \v] .rtj d 3 )k 3 



Yi = 



de d 
' dd. 



where 



Y 2 
Y 3 



1 



? 



fciMi 

k 2 [u]l 



(3) 



= ^3 ([u] 3 )+ , 




Figure 1: The mesomodel entities 



where e^ is the strain energy of the interface per unit area and [u]i is the i th component 
of field [u]. We assume that the damage variables are functions of a single quantity: 
the maximum Y\ t over time of a combination of the energy release rates lii r , r < t: 
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Thus, the evolution law is defined by: 



d\ = d>2 = ds = w(Y") 



where, in general, k;(Y") 



(4) 



(5) 



1 



This model has the advantage, using a single damage variable, to be able to recover 
different critical energy release rates for pure mode loading as well as to fit to classical 
mixed mode criteria. For example, setting a = 2 leads to the classical quadratic energy 
criteria. 



3 Main features of the method 

As presented in Figure [2| the structure is split into substructures which are connected 
by interfaces that have mechanical behaviors. In a second stage, the substructures are 
gathered into "super-substructures" which are assigned to independent processors. In 
order to ensure the scalability of the method, the interfaces are the support of a macro 
(homogenized) problem, the resolution of which requires the solution of a third scale 
(supermacro) problem on super-interfaces. 



3.1 Micro-macro scale separation 

The iterative LaTIn algorithm, which is designed to non-linear problems, is here ap- 
plied to the solution of the geometrically non-linear substructured reference problem 
with "material" non-linearities (damage, contact) localized on the interfaces. The finite 
element method is used to discretize the governing equations. The aim of the method 
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Figure 2: Decomposition of the problem and associated three computational scales 



is to find the subdomains fields u Eq (displacement), n (Piola-Kirchhoff stress), and 
the interface fields W r (displacement), F p (interforces), where index E ranges over 

^ ^ 

all substructures and index means that the quantity is given in the reference config- 
uration (total Lagrangian approach). 

In order to ensure the scalability of the method, a global and linear coarse grid 
problem is solved. The definition of the macroscopic fields required to construct this 
problem is done on the interface unknowns only. Whichever the interface behavior, 
action reaction principle F_ E + F_ E , = is verified; the principle of the macro-problem 
is to enforce this equation partially at any time of the iterative solution process: 
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where displacement macrospace W E is a parameter of the method as well as its dual 
T E . The definition of the macro-spaces is done through a projector II, which is chosen 
L 2 -orthogonal (so that the same projector is used by the traction and displacement 
fields). We then have the following splitting of the interface quantities: 
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(7) 
Numerical tests showed that in order to ensure the numerical scalability of the 
method the macro-basis should extract at least the linear part of the interface forces. 
Indeed, this macro-space contains the part of the interface fields with the highest wave- 
length. Consequently, according to the Saint-Venant principle, the micro complement 
only has a local influence. 

3.2 Iterative algorithm 

It is now possible to split all the equations of the system into two groups: 

• non-linear equations in the substructures and macroscopic admissibility of inter- 
faces, whose solutions are elements of the manifold A^: 

- non-linear kinematic admissibility of the substructures 

K Eo = \ (^ u Eo + %u Eo + % u Eo % u Eo ) , on n Eo (8) 

UE omEo =W EolrEo ,onT EoE , (9) 
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non-linear static admissibility of the substructures 
V(u Eo ,W* Eo )eU Eo xW Eo such that u EolmEo =Wl o]TEo , 

k Eo ■ k(u Eo ) dn 

= / PE l d -V*E d£lo+ I F Eg -W Eo dT (10) 

- behavior of the substructures 

dip 
=e = qe~ > on n *o , (11) 

- macroscopic admissibility of the interfaces (|6| . 

• local (non-linear) equations in the interfaces whose solutions belong to the man- 
ifold T: 

- interface behavior (perfect, contact and cohesive interfaces) and boundary 
conditions (see Section 2). 



The solution s re f 



{u Eo ,^ W F 



is such that: 
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Note that in the case of small perturbations, the admissibility equations are linear. In 
this case, the manifold Ad is an affine space. 

The resolution scheme consists in seeking the solution s re f alternatively in these 
two manifolds: first, one finds a solution s n in Ad, then a solution s n+ i in T. In 
order for the two problems to be well-posed, search directions E + and E~ linking the 
solutions s and s through the iterative process are introduced (see Figure [3| . 

• At the global stage, starting from fields (W Eo , F Eq ) satisfying the interface equa- 
tions (manifold T), we seek fields (W E ,F_ E ) in Ad (satisfying the subdomain 
equations) using the following closing relation: 



k Eo (W Eo 



w Eo ) + ( 



F T 



J 



F Eo \=0 



(13) 



k E is the search direction associated to the global stage. 

At the local stage, starting from fields (W w iKt?) satisfying the subdomain 

equations (belonging to manifold Ad), we seek fields (W E ,F E ) satisfying the 
interface equations (manifold T), using the following closing relation: 
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k E is the search direction associated to the local stage. 
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Figure 3: Schematic representation of the iterative LaTIn algorithm 
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Figure 4: LaTIn convergence curves for different numbers of macro iterations 



3.3 Third scale 

The macro-problem being linear, discrete and sparse, it can be efficiently solved by a 
balancing domain decomposition algorithm [41J. To do so, substructures are gathered 
into super-substructures. The macro-problem is condensed at the super-interfaces and 
this condensed problem is solved by a conjugate gradient algorithm. The balance of 
super-structures requires the global solution of a third scale problem whose size is the 
number of rigid body motions of floating super-substructures. 

Figure Q shows the convergence rate of the LaTIn algorithm when the conjugate 
gradient scheme for the condensed macroproblem is stopped after a fixed number of 
iterations. The test case is the holed plate under traction loading with the super- 
substructuring pattern given in Figure (|5|. This problem is 3.4 10 6 degrees of freedom 
large, it is distributed in 520 subdomains with 1,350 interfaces leading to a 12 10 3 - 
unknown macro-problem. 8 super-substructures are employed so that the condensed 
macro-problem on the super-interfaces has 1, 746 unknowns and the dimension of the 
third scale (associated to the rigid body motions of the super-substructures) is only 
36. 

It appears clearly that only very few iterations of the conjugate gradient scheme 
are required to get the necessary part of the macro-displacement Lagrange multiplier 
leading to the multiscale convergence rate of the LaTIn algorithm. Typically, the 
algorithm is stopped when the residual error (normalized by the initial error) falls 
below 10 _1 . The admissibility of the macro- forces is thus enforced on a third level, 
which is sufficient to determine the part of the solution which needs to be transmitted 
at each iteration of the resolution. In that case the time required to solve the macro- 
problem is divided by a factor 100. 

Figure [6] presents a larger case where a direct computation of the macroproblem 




Figure 5: Different levels of substructuring 



would be too expensive. This 12 10 6 dof problem is distributed on 10, 960 substructures 
and 32 10 3 interfaces, leading to a 300 10 3 macro d.o.f, the super-substructuring leads 
to a 40 10 3 condensed dof problems and 150 super-macro d.o.f. Figure [7| shows the 
extension of delamination in the final converged state. 




Figure 6: Bolted joint, substructuring and super-substructuring 

The previous examples allow us to make an important remark related to our vision 
of the problem of large-scale simulation of delamination. The reader will notice that 
we use a homogeneously fine finite element discretization to approximate the displace- 
ment field in the structure. Specifically, the mesh is sufficiently fine to handle the 
propagation of a delamination front at an arbitrary location in the structure, which 
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Figure 7: Final delamination state of the bolted joint 




Figure 8: Relocalization box around the crack tip 

permits to handle multiple and arbitrary large delamination fronts without remeshing. 
The fact that we are not willing to remesh is a choice driven by the observation that 
error estimation, mesh adapt ivity and field transfer are, for parallel problems and for 
problems with softening behaviors, a subject of research in themselves and no ma- 
ture tool has emerged from pioneer investigations which can be found, for instance, 
in [331 ESI El| EJ [43]. Of course, the price that we pay for the expected robustness 
associated to this choice is a large number of arguably unnecessary degrees of freedom. 

4 Relocalization for the propagation of delamination 

Let us consider a simple DCB case, as illustrated in Figure [8J Figure [9] shows that 
the convergence of the classical approach (blue bars) is strongly deteriorated when 
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Time step 



Figure 9: Performance of the relocalization technique w.r.t. the size of 
the relocalization box 

delamination starts to propagate (time step 3). Indeed, in average 10 times more 
iterations (and CPU time) are required. This drop in the convergence rate occurring 
when the cracks propagates can be explained by two driving factors: 

• the singularity near the tip of the crack is very poorly represented by linear macro 
quantities. Therefore, the complementary part is not localized enough and many 
iterations are required to transmit them through the structure; 

• the strong non-linearities require many iterations to converge, in particular the 
prediction of the location of the crack tip at a given time step requires the quasi- 
convergence of a large number of consecutive equilibrium states. 

Using non-linear relocalization in a box containing the crack tip and the process zone 
is a way to deal with those two problems: after each linear stage, a subproblem is 
extracted around the crack tip (within the so-called relocalization box whose size is 
a parameter) as shown in Figure [81 This smaller problem is solved non-linearly with 
Robin boundary conditions satisfying the macro-equilibrium of the structure. 

Being far enough from the singularity, the macro quantities at the boundary of the 
box represent the global fields correctly and are then relevant for the prediction of the 
non-linear evolution of the crack. The resulting modifications only have a small-scale 
effect which can be transmitted by the upcoming local stage. 

The other bars in Figure [9] show the influence of the size of the relocalization box on 
the convergence. Of course the larger the box, the best the filtering of the singularities 
on the boundary of the box, and the fastest the convergence. The counterpart is that 
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the larger the box, the more expensive the non-linear computations associated to the 
relocalized solutions, which are obtained by a LaTIn method restricted to the box, are. 

5 Further improvements for delamination/buckling and 
contact 

5.1 Taking into account the structure's slenderness 

Optimal search direction of one subdomain (for the global step) is known to be the 
Schur complement of the rest of the structure. In general a very coarse approximation, 
in practice a scalar, is sufficient to ensure good convergence of the LaTIn method. In 
[47J it was shown that for plates, the slenderness induces a structural anisotropy that 
has to be taken into account by the search directions. A dimensional analysis leads 
to the following relationship between the longitudinal (n) and transverse (t) search 
directions (L is the characteristic length of the structure, ho is its thickness): 



m 



m 



h 



Here superscript m stands for the micro part of the search direction. The macro part 
of the search direction has to be set to a large value on perfect interfaces so that the 
macro continuity of displacements is almost satisfied throughout the iterative process. 
Figure [To] shows that these adaptations of the search directions (i.e. anisotropic 
search directions) are necessary to ensure the scalability of the method. 
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Figure 10: Influence of the search directions and of the substructuration in a bending test 
case 
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5.2 Update of the search directions 

In the context of buckling and contact, the state of the structure can have large in- 
fluence on its response. This implies that search directions which are optimal in one 
configuration can be unadapted to others, leading to a drop in the convergence rate. 
Strategies to update the search directions according to the status of the interfaces must 
then be settled. Typically for cohesive interfaces, the optimal value is related to the 
current stiffness of the interface (which depends on the damage state) ; for fully delam- 
inated interfaces, search directions should be large when contact occurs whereas they 
should be null when the crack is open. 

Figure [TT| presents the effect of the updating of the search direction according to the 
damage state. In this simple DCB case, all strategies converge. Of course, the more 
frequent the updates the less iterations, but each update requires recomputation and 
factorization of the operators so that a compromise must be found to achieve minimal 
CPU time. 

Figure [12] shows the evolution of the error for one time step of a opening contact 
problem: the plate is initially delaminated in its middle (large closed crack) and a 
compression loading is imposed. When the optimal search direction is chosen a priori, 
only 10 iterations are required. When the wrong one is chosen, the algorithm converges 
to a non-physical solution. Starting from an incorrect initial guess, a simple update 
procedure (every 10 iterations the parameters are adapted to current contact state) 
enables to converge to the solution. Note that the convergence is slow because of the 
unusual extension of the initial contact zone. 

5.3 Pre-delaminated plate sensitive to imperfection 

Let us consider the example of a pre-delaminated plate submitted to compression. The 
example is designed to be sensitive to imperfection. Indeed two buckling modes can be 



excited (Figure 13): a global buckling mode (closed interface - mode II delamination) 
and a local buckling (open interface - mode I delamination). The problem is split 
into 1280 substructures leading to 3248 interfaces with about 1.2 millions of d.o.f. The 
macro-problem involves 30,000 d.o.f. when the super-macro one only involves 132 d.o.f. 
When the imperfection is small, the response of the buckling mode of the plate is 
a global one (mode II delamination) and the delamination starts to propagate for a 



transversal displacement of about 0.75 mm (see Figure 14 - mode II curve). 

When the level of the perturbation is at least ten times greater {F_ d = 2 N), a local 
buckling mode is excited, which leads to the propagation of the delamination for a 



transversal displacement of about 0.36 mm (see Figure 14 - mode I curve). Figure 15 



shows the configuration and the crack propagation of the last time step. 

6 Example of multiple delamination of a laminated 
plate in compression 

The following example has been chosen because it involves all the difficulties discussed 
in this last section. Here, the delamination between plies and the possibility of contact 
and separation between delaminated surfaces in large displacements are studied (others 
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Figure 11: Iterations and CPU time to convergence when using different updating strategies 



possibles scenarios like transverse cracking are obviously ignored by the model). It is 
a [0°/90°] s composite plate, submitted to compression (prescribed displacement u d ) 
and pre-delaminated (with the same initial length) on the central interface and on 
the upper one. A normal perturbation is applied F d at the level of the top ply (see 



Figure 16). The layers are assumed to be elastic and orthotropic with the following 
characteristics in the referenced basis (1- fiber direction, 2- transverse direction in the 
ply, 3- normal direction): E 1 = 185, 500 MPa, E 2 = E 3 = 9, 900 MPa, u 12 = viz = 0.34, 
za>3 = 0.5, G12 = G13 = 6, 160 MPa, G 23 = 3, 080 MPa. The problem is split into 1,280 
substructures leading to 3,248 interfaces, 2 millions of d.o.f. are involved. The macro- 
problem involves 30,000 d.o.f. when the super-macro one only involves 168 d.o.f. The 
time interval is discretized into 120 times steps and 30 processors have been used. On 
Figure RTFl the force-normal displacement responses for the three groups of layers in 



the middle of the plate are plotted and Figure [18] shows three typical configurations. 
The top ply buckles first for an applied load of 80 N leading to a non symmetric 
configuration (configuration A). In this configuration a loss of contact between the 
two upper plies and the lower one is observed. When the applied load reaches 100 
N, the central and lower plies buckle which leads to a recontact between the upper 
plies (configuration B). Finally, the configuration C is symmetric. For a value of the 
transverse displacement of about 0.5 mm, the delamination starts to propagate in the 
central interface leading to a decrease of the reaction forces. 
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Figure 12: Evolution of the error when using different updating strategies 

7 Conclusion 

In order to circumvent the numerous problems associated with the large scale simula- 
tion of delamination problems, we have developed a multi-scale framework, wherein a 
Domain Decomposition approach is optimized to deal with the severe non-linearities. 
A lot of work is still required to make the approach more generic. The problem of code 
development is also an issue that needs to be taken into account when developing such 
dedicated methods. The challenge we would like to address is the one of the Virtual 
Delamination Testing, in a fully integrated sense. Nevertheless, those studies have been 
conducted by restricting the potential damage modes of the laminate to delamination. 
However, it is known that in many cases, like in the case of low- velocity impacts or in 
the case of the analysis of the tolerance to defects, lumping the damages onto the inter- 
faces cannot reproduce the delamination patterns adequately. In fact, these patterns 
are often driven by the interaction between transverse cracking and delamination. 

Two main approaches are followed in this more general case. In [40J a micro analysis 
leading to a meso-scale model where a non-local coupling between the interface and 
the adjacent plies associated with transverse-cracking / delamination interaction has 
been developed. A validation of this approach by making use of open-holed tests, 
which exhibits very different delamination patterns due to scale effect [23J can be 
found in [lj. The non- locality of the ply model requires to have all the information 
throughout the thickness of a given ply within a unique sub-domain. This is not a 
problem with the physical decomposition at the meso-scale proposed in the paper. 
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Figure 13: Pre-delaminated plate with associated potential buckling modes 
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Figure 14: Load displacement curves depending on the imperfection 



However, it becomes an issue if one makes use of an automatic decomposition of the 
structure, as was addressed in [8J. Another approach aims at introducing the transverse 
cracks explicitly, or at least the transverse cracks that have a strong influence on 
the delamination pattern [50J by means of X-FEM like technique. In the framework 
proposed in this paper, following this direction requires to development an hybrid 
multi-scale X-FEM method, as was proposed in |24|. However, the complexity and 
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Figure 15: Final state of the plate in case of local buckling 

initially cracked interfaces (contact) 
cohesive interface 




top layer 
upper 90 layer 

J lower layers 



central section 



Figure 16: Configuration of the specimen 
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Figure 17: Displacement of three sets of plies during the experiment 
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Figure 18: Three configurations. A: all plies separated, B: recontact, C: ruin 



highly- intrusive character of such an approach is an issue in terms of code development. 
A promising alternative is to perform a local enrichment by means of a non-intrusive 
coupling of the structural solver with specialized external solvers dedicated to specific 
failure modes, as was proposed in [21]. For instance, one can couple the tools presented 
in this paper with the one developed to simulate a micro-model where all the cracks 
are introduced explicitly |51j . This idea was successfully applied in [14J. 
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